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Abstract We present a radiation hydrodynamics simulation of the formation of a 
massive star using a Monte Carlo treatment for the radiation field. We find that 
strong, high speed bipolar cavities are driven by the radiation from the protostar, and 
that accretion occurs stochastically from a circumstellar disc. We have computed 
spectral energy distributions and images at each timestep, which may in future be 
used to compare our models with photometric, spectroscopic, and interferometric 
observations of young massive stellar objects. 



1 Introduction 

The deposition of substantial momentum from the radiation field into the accreting 
material of the protostellar envelope creates a serious obstacle for theories of mas- 
sive star formation that adopt spherical accretion e.g. [1 1. Extending the models to 
multiple dimensions and the inclusion of rotation leads to a scenario in which the 
protostar can accrete via a disc, which subtends a smaller solid angle at the proto- 
stellar surface and thus intercepts a smaller fraction of the radiation e.g. ID. 

Recent work by Krumholz et al. used adaptive mesh refinement hydrodynamics 
coupled with a flux limited diffusion (FLD) approximation for the radiative transfer 
(RT) in order to simulate the formation of a massive binary 1 3 1 . They found accretion 
occurred via a radiatively-driven Rayleigh Taylor instability that allowed material to 
penetrate the radiation-driven cavities above and below the protostars. Simulations 
by Kuiper et al. ||4| using a hybrid RT method 15| have cast doubt on this scenario 
- they claim that the grey FLD approximation leads to a substantial underestimate 
of the driving force of the radiation field (see also Kuiper's contribution to these 
proceedings). 
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Clearly the level of detail with which the radiation transport is followed is an 
important consideration when simulating radiation feedback. We have therefore de- 
veloped our own radiation-hydrodynamics code |6| which is based on the TORUS 
Monte Carlo RT code QISIIH. The advantage of this approach is that RT is fully 
polychromatic and includes both dust and atomic microphysics, and although the 
method is unavoidably computationally intensive, the overhead is reduced to ac- 
ceptable levels by the use of aggressive paralleUzation. 



2 The Monte Carlo method and its implementation 

The TORUS code is based on the radiative-equilibrium method of ifTOl in which the 
stellar luminosity is divided up into a large number of photon packets, which prop- 
agate through the computational domain undergoing absorptions and scatterings. 
When a photon packet is absorbed it is immediately re-emitted with a frequency 
found from random sampling of the emissivity spectrum at that point in the grid. By 
tracking the paths of the packets through a grid cell one can obtain a Monte Carlo 
estimate for the energy density in that cell, and hence the mean intensity of the 
radiation field, which in turn enables the computation of absorption rates, photoion- 
ization rates etc. The radiation force on each cell can be computed using similar 
methods. 

The primary drawback in adopting such a detailed treatment of the radiation field 
is the computational effort involved. Fortunately the Monte Carlo method, in which 
the photon packets are essentially independent events, may be straightforwardly and 
effectively parallelized. The top level of parallelization involves domain decompos- 
ing the octree that stores the adaptive mesh. The main bottleneck is the communica- 
tion overhead when passing photon packets between threads, and we optimize this 
by passing stacks of photon packet data between the threads rather than individual 
packets. A further level of parallelization is achieved by having many sets of iden- 
tical computational domains, over which the photon packet loop is split, with the 
results derived from the radiation calculation (radiation momenta, absorption rate 
estimators etc) collated and returned to all the sets at the end of each iteration. Us- 
ing these parallelization methods it is possible to reduce the time for the radiation 
field calculation down to a level comparable with that of the hydrodynamics step. 



3 A test calculation 

We adopted the initial conditions from one of the models presented by f4\ in order 
to perform a demonstration calculation. This model consists of a 100 M0, 0.1 pc 
radius cloud with an r^'^ density profile and a small amount (£2 = 5 x 10^'^ s^^) 
of solid body rotation. The effective temperature of the protostar and its luminosity 
were determined by interpolating by mass in the protostellar evolution models of 
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Hosokawa & Omukai ifTTI . The photospheric spectrum was composed of an appro- 
priate Kurucz model plus a blackbody accretion spectrum. A 128^ regular mesh was 
used for this test calculation. 

The final state of the model (after 50 000 years) is plotted in Fig.[T] The radiation- 
driven bipolar cavities are clearly seen, and the outflow has reached the boundary of 
the computational domain. The development and speed of the cavities is similar to 
that found by Kuiper et al., but the morphology is different, principally due to the 
resolution. We have also computed SEDs at each timestep, and plot these in Fig. |2] 
Note that 60° and 90° inclination SEDs show typical Class I characteristics, and 
the attenuated protostellar SED is observed only directly along the outflow cavities. 
The evolution of the accretion onto the protostar can be divided into two distinct 
phases: At early times the accretion is smooth, and from the envelope, but after 
^-^20000 years the luminosity of the protostar overtakes the accretion luminosity 
(this is dictated by the evolutionary tracks) and the radiation-pressure starts to drive 
the polar outflow cavities, at which point accretion is occurring stochastically via 
the protostellar disc. (Note that this two stage accretion evolution is also present 
in the hybrid RT calculations of Kuiper et al.) At the end of the calculation the 
protostar has a mass of 28 Mq and an effective temperature of 30 000 K, at which 
point ionization and stellar wind feedback are potentially important. Although a 
small ionized zone is formed close to the protostar towards the end of the calculation 
(with properties that would identify it as an ultra-compact H II region) we note that 
we have yet to implement a self-consistent method for treating the dust sublimation 
and advection. 
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Fig. 1 A slice through the 
rotation axis of the test sim- 
ulation (spatial units are lO'" 
cm). The greyscale shows 
logarithmic density scaled be- 
tween 10^" gcm^' (white) 
and 10^'^ gcm^' (black). 
The arrows denote velocity, 
with the longest arrows cor- 
responding to speeds of 100 
kms^'. 
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Fig. 2 The results of the test calculation. Spectral energy distributions (ergs^' cm^^ at a distance 
of 1 kpc) at inclinations of 0, 60 and 90 degrees are shown in the top-left figure. The mass of 
the protostar as a function of time is plotted in the top-right figure. The effective temperature and 
luminosity of the protostar is shown in the bottom left figure (note the temperature increases to the 
left). The bottom-right figure shows the accretion rate onto the protostar as a function of time. 
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